The Epidermal Transcriptome Analysis of a Novel c.639_642dup LORICRIN Variant-Delineation of the Loricrin Keratoderma Pathology

Loricrin keratoderma (LK) is a rare autosomal dominant genodermatosis caused by LORICRIN gene mutations. The pathogenesis of the disease is not yet fully understood. So far, only 10 pathogenic variants in LORICRIN have been described, with all of them but one being deletions or insertions. The significance of rare nonsense variants remains unclear. Furthermore, no data regarding the RNA expression in affected patients are available. The aim of this study is to describe the two variants in the LORICRIN gene found in two distinct families: the novel pathogenic variant c.639_642dup and a rare c.10C > T (p.Gln4Ter) of unknown significance. We also present the results of the transcriptome analysis of the lesional loricrin keratoderma epidermis of a patient with c.639_642dup. We show that in the LK lesion, the genes associated with epidermis development and keratocyte differentiation are upregulated, while genes engaged in cell adhesion, differentiation developmental processes, ion homeostasis and transport, signaling and cell communication are downregulated. In the context of the p.Gln4Ter clinical significance evaluation, we provide data indicating that LORICRIN haploinsufficiency has no skin consequences. Our results give further insight into the pathogenesis of LK, which may have therapeutic implications in the future and important significance in the context of genetic counseling.


Background
Loricrin keratoderma (LK, MIM 604117, Vohwinkel syndrome with ichthyosis, VS) is a rare autosomal dominant genodermatosis caused by pathogenic variants in the LORICRIN gene. LORICRIN encodes one of the key proteins conferring the insolubility, mechanical resistance and impermeability of the epidermal barrier [1] Hydrophobic and insoluble loricrin is expressed in the orthokeratinizing epithelia, except for internal ones. In the skin, its synthesis occurs in the upper layer of epidermisstratum granulosum (SG). Loricrin is involved in cytoskeleton stabilization forming crosslinks within and between the proteins, and in the formation of the cornified cell envelope (CE), being the most abundant (commonly > 70%) protein there [2,3].
Consequently, the clinical symptoms of LK patients are related to the skin surface and comprise the following: ichthyosis; palmoplantar keratoderma (PPK), often with honeycomb pattern, pseudoainhum and/or amputation; knuckle pads; and collodion membrane at birth [4]. Of note, the clinical symptoms are heterogenous and may differ even among relatives [5].
The data concerning LK are limited. Only 10 pathogenic variants in 21 affected families (overall 106 patients) were described in the literature so far [4]. Moreover, all of them, apart from one substitution, are deletions/insertions. The only pathogenic missense variant known so far was identified as causative in late-onset loricrin keratoderma [6]. The clinical significance of the other variant types remains questionable.
The aim of the study is to describe two variants in the LORICRIN gene that were found during diagnostic procedures of cornification disorders. The first one, c.639_642dup (p.Thr215GlyfsTer122), is a novel pathogenic variant detected in a family with autosomal dominant hyperkeratosis, for which we also present the results of a transcriptomic analysis. This is the first transcriptomic analysis of a loricrin keratoderma lesion. Another variant, the rare c.10C > T (p.Gln4Ter), was detected in the other family as a secondary finding of unknown significance. Considering the highly limited data on the clinical significance of LORICRIN premature stop codon (PTC) variants, we provide data showing that p.Gln4Ter leading to a premature stop codon has no skin consequences. This has important significance in the context of genetic counseling.

Experimental Design
All patients gave informed consent to participate in the study.

Patients
Family 1: The family (two daughters and their father, Figure 1A) was referred to genetic counselling because of hyperkeratosis of the palms and soles and a clinical diagnosis of ichthyosis. The clinical symptoms were manifested by ichthyosiform dermatosis, diffuse generalized. In the girls, the symptoms were noted at birth. Then, palmoplantar keratoderma occurred at the age of 2-3 months. The course of the disease varied, with occasional exacerbations. The improvement was noted after the use of emollients. Occasionally, the transgradient extension of hyperkeratosis onto the wrists and on the bends of the elbows and knees was present, pseudoainhum was not observed and, according to the patient, was also absent in the other affected family members. The honeycomb pattern of PPK was negative during clinical evaluation and, according to the mother, had not been observed before. The keratoderma was neither painful nor inflammatory.
Family 2: The proband was a girl born from an uneventful pregnancy at 38 weeks of gestation (birth weight 2880 g, Apgar score: 7). She had clinical recognition of autosomal recessive congenital ichthyosis (ARCI), due to a homozygous pathogenic variant c.1562A>G (p.Tyr521Cys) in ALOX12B. The symptoms of ARCI were typical (collodion baby; laterin-life dryness of the face skin and, less intensive, of the whole body; fingers and toes contractures; erythema; stiff and cracking skin of the hands and feet; slight psoriasis lesions on the knees and elbows), no nail and hair disturbances were present and the teeth appeared normally, though a slight yellow discoloration of permanent teeth was observed ( Figure 2).
The LORICRIN variant p.Gln4Ter was detected as a secondary finding during a molecular test. A segregation analysis has shown that the variant was inherited from the patient's father. A dermatological evaluation of the father did not reveal any skin symptoms at the age of 41; only dystrophic nails were present and massive caries (currently with upper teeth dentures) from the age of 20. Family 2: The proband was a girl born from an uneventful pregnancy at 38 weeks of gestation (birth weight 2880g, Apgar score: 7). She had clinical recognition of autosomal recessive congenital ichthyosis (ARCI), due to a homozygous pathogenic variant c.1562A>G (p.Tyr521Cys) in ALOX12B. The symptoms of ARCI were typical (collodion baby; later-in-life dryness of the face skin and, less intensive, of the whole body; fingers and toes contractures; erythema; stiff and cracking skin of the hands and feet; slight psoriasis lesions on the knees and elbows), no nail and hair disturbances were present and the teeth appeared normally, though a slight yellow discoloration of permanent teeth was observed ( Figure 2).   Family 2: The proband was a girl born from an uneventful pregnancy at 38 weeks of gestation (birth weight 2880g, Apgar score: 7). She had clinical recognition of autosomal recessive congenital ichthyosis (ARCI), due to a homozygous pathogenic variant c.1562A>G (p.Tyr521Cys) in ALOX12B. The symptoms of ARCI were typical (collodion baby; later-in-life dryness of the face skin and, less intensive, of the whole body; fingers and toes contractures; erythema; stiff and cracking skin of the hands and feet; slight psoriasis lesions on the knees and elbows), no nail and hair disturbances were present and the teeth appeared normally, though a slight yellow discoloration of permanent teeth was observed ( Figure 2).

Genetic Analysis
We identified a novel variant in the LORICRIN gene: c.639_642dup (p.Thr215GlyfsTer122) in family 1 ( Figure 1A). The pathogenicity status was scored as likely pathogenic (LP) according to the American College of Medical Genetics (ACMG) classification. Importantly, similarly to other LORICRIN pathogenic variants reported so far, the c.639_642dup caused delayed translation termination and introduced an arginine and leucine reach sequence. Thus, the diagnosis of loricrin keratoderma was established.

Transcriptome Analysis of the Probant vs. Control
The 15,210 genes with distinct ensemble identification (ID) and more than five counts in each sample were detected. Considering the fact that the data analysis was largely limited and included a single patient vs. single control analysis, we highly strengthened the differentially expressed genes (DEG) parameters to the absolute value of logarithm fold change (|logFC|) > 3 and logarithm of counts per million reads (logCPM) > 1, resulting in 1722 genes. Among them, 276 genes were upregulated (logFC between 3.05 and 12.7) and 1445 downregulated (logFC between −3.0 and −14.45). However, only in 10 and 53 genes, respectively, was the statistical difference significant (p-value < 0.005) ( Table 1). With respect to ontology, genes-encoding proteins involved in epidermis development and keratinocyte differentiation were mainly upregulated. In turn, those engaged in, i.e., cell adhesion, developmental processes and anatomical structure morphogenesis, cellular ion homeostasis and transport, cell differentiation, regulation of signaling and cell communication were downregulated ( Table 2).

Discussion
In this study, we described the novel LORICRIN gene pathogenic variant: c.639_642dup with the first transcriptome analysis of lesional loricrin keratoderma epidermis, and a rare p.Gln4Ter variant in the same gene, as evidence that the haploinsufficiency of loricrin does not cause skin symptoms of LK.
The in silico prediction showed that the consequence of c.639_642dup on the protein level is a generation of a sequence rich in basic amino acids, mainly arginine. It has already been proven that all the other known frameshifts in the C-part of the loricrin also lead to the formation of arginine-rich regions generating nuclear localization signals (NLS) [7]. Indeed, such loricrin derivative mutated proteins were found to be deposited in the nucleus and distort epidermal differentiation [5,8]. This was also observed in a mouse model of LK, where mutated loricrin was almost exclusively present in the nucleus. This, in fact, was further proven to be an LK-causative factor. It was also shown that the LK phenotype of transgenic mice was more severe in the absence of wild type loricrin [8].
Next-generation sequencing technologies (NGS) enabled robust progress in the genetics of the disorders of cornification. While DNA sequencing has already revealed a plethora of disease-causing variants, showing great heterogeneity in the molecular basis of these diseases, RNA sequencing data from these patients are rather limited. Nevertheless, a few studies have already shown that transcriptome analyses may be crucial for obtaining deeper insight into the pathophysiology of the cornification disorders. However, according to our knowledge, no data on the gene expression in loricrin keratoderma patients are available, probably due to the rarity of this disorder.
Herein, we showed the results of the transcriptome analysis performed using mRNA isolated from the lesion epidermis of the patient with heterozygous novel variant c.639_642dup. In total, 1722 genes were differentially expressed, of which 276 genes were upregulated and 1445 downregulated. However, only 10 and 53 genes reached statistical significance, respectively.
The HRNR-encoding hornein was the most upregulated gene. This gene is located on chromosome 1q21 within the human epidermal differentiation complex (EDC). Hornein belongs to S100 fused-type proteins (SFTPs) and is involved in the cornified epithelium formation [9]. Furthermore, hornein has an antimicrobial activity as the source of cationic intrinsically disordered antimicrobial peptides (CIDAMPs) [10,11]. It has previously been shown that HRNR mRNA expression increased transiently in cultured human epidermal keratinocytes during Ca 2+ -dependent differentiation [12]. Of note, Rice et al. and Kim et al. have shown that in healthy people, the HRNR is preferably expressed in palmoplantar skin compared to other regions [13,14].
So far, the HRNR gene was mainly analyzed in the context of the other skin diseases: psoriasis and atopic dermatitis (AD), where barrier defects occur as well, but due to distinct immunogenetic factors. The HRNR transcripts were detected in regenerating human skin after wounding in the periphery regions of psoriatic lesions [15]. Moreover, the hornein immunoreactivity in the lesions, but not in the healthy skin, of psoriasis and atopic dermatitis patients was also diminished in another study [12]. Furthermore, Henry et al. showed that the expression was lower also in the healthy skin of AD patients. The authors demonstrated that hornein is a component in a cornified envelope (CE) and suggested that it plays a role in the alterations in the CE and in the abnormality of the AD epidermal barrier [16].
Just recently, Makino et al. checked the HRNR expression by immunostaining in skin lesions from patients affected by hyperkeratosis-associated diseases (ichthyosis vulgaris, epidermolytic ichthyosis (EI), Darier's disease, lichen planus, pustulosis et plantaris, actinic keratosis and seborrheic keratosis). The increased expression was detected in lichen planus and pustulosis et plantaris, followed by an irregular signal pattern in epidermolytic ichthyosis and actinic keratosis. In the remaining diseases (ichthyosis vulgaris, Darier's disease and seborrheic keratosis), the expression was decreased. Thus, in light of our results and those mentioned above, further studies are needed to evaluate the hornein involvement in epidermal barrier restoration [17].
Among the other top 10 upregulated genes, we detected a few more associated with barrier formation: LCE3D (late cornified envelope protein 3d), KRT9 (keratin 9) and CDSN (corneodesmosin). Those genes were also found to be upregulated in the other types of ichthyoses [18]. Specifically, LCE3D was also found to be upregulated in the other diseases with keratoderma: Pachyonychia Congenita and Curth-Macklin ichthyosis [19,20].
Due to the fact that our analysis consisted of only one patient and one control, the statistical analysis was very limited. Therefore, we also focused on the genes that had logFC over 3.0 or below −3.0 and logCPM > 1, irrespective of the p-value. In this group, several others had induced expression as well, including the IL-17/TNF-α-associated molecules IL36G and S100A9. Of note, previous studies also showed that the Th17 pathway is induced in various forms of ichthyosis, which proves that in terms of immune response, ichthyoses resemble psoriasis. Hence, novel therapies using IL-17 may be deliberated in the future [18,21]. Overall, among the 276 upregulated genes, those associated with epithelium development, keratocyte differentiation and keratinization were the most represented.
It has been shown that apart from some commonly expressed genes, different ichthyoses vary in terms of gene expression. This is reflected even by the numbers of DEG (patient's lesions vs. control) in different disorders. Malik et al. showed that in the Netherton syndrome patient, the number was relatively low: 63 upregulated and 33 downregulated DEGs comparing to epidermolytic ichthyosis (EI), where the number of DEGs was 223 and 150, respectively. Furthermore, Kim et al. identified lipid metabolism and barrier junction genes to be downregulated in four common ichthyosis types, which were less pronounced in EI [21]. Furthermore, Malik et al. proved that the expression of lipid metabolism genes was diminished in lamellar ichthyosis (LI) patients, but not, or to a lesser extent, in EI [18]. This phenomenon may result from the distinct molecular basis of those disorders: LI is mainly caused by mutations in genes involved in lipid metabolism, while in EI, pathogenic variants in structural keratins 1 and 10 are causative. Finally, when we compared the genes downregulated in our patient with those published by Ortega-Recalde et al. in Curth-Macklin ichthyosis, two were concordantly downregulated (PHYHIP, PAMR1), while DCD, FABP4, PLIN1, SCGB1D2, SCGB2A2, ADIPOQ, G0S2, KRT19 and MUCL1, downregulated in our case, were upregulated in Curth-Macklin ichthyosis.
Among the mostly downregulated 53 genes, we found a few involved in lipid metabolism, e.g., PLAAT3, FABP4, PLIN1/4 and PLNPLA6. However, a gene ontology analysis performed in the wider context showed that the majority of 1445 genes detected by us are involved in cell adhesion developmental processes and anatomical structure morphogenesis, cellular ion homeostasis and transport, cell differentiation, regulation of signaling and cell communication. This finding is in line with the molecular pathomechanism of loricrin keratoderma, which, as already mentioned before, comprises the nuclear deposition of mutated loricrin and the dysregulation of keratinocyte differentiation. Of note, once we compared the 53 mostly downregulated genes of our patient with the DEG profile of atopic dermatitis and psoriasis presented by Malik et al., 14 also had a diminished expression in AD and 15 in psoriasis, whereas only a few (1)(2)(3)(4)(5) overlapped with other ichthyoses [18].
Since the results, to our knowledge, are the first transcriptome analysis of LK lesion and were performed on the one patient only, replicative studies are needed. Nevertheless, our results provide novel insight into the pathogenesis of the disease and may have therapeutic implications in the future.
Another issue raised by us concerns the clinical significance of nonsense variants in LORICRIN. It has been shown that transgenic mice with one copy of the loricrin gene are phenotypically normal [22]. However, as far as we can tell, there are no phenotypic descriptions of humans with one functional copy of the LORICRIN gene available thus far. Hence, the genetic counseling in such cases may be ambiguous.
In the ClinVar database, one premature stop codon (PTC) variant [NM_000427.3:c.624C > G (p.Tyr208Ter), ID: 1324671] is recorded and is assigned as likely pathogenic. On the contrary, in the SNP database (SNPdb), 13 nonsense variants are recorded, with the frequency ranging from 0 to 0.00004, according to the GnomAD or Kaviar databases. None of the variants were detected in homozygosity and each was classified as VUS according to the ACMG [23] classification. The variant c.10C > T (p.Gln4Ter) detected in family 2 is also recorded in SNPdb and was found in 2 out of 231 412 GnomAD alleles of European, non-Finnish ancestry.
The proband of family 2 was diagnosed as autosomal recessive congenital ichthyosis ARCI with ALOX12B biallelic mutations; therefore, it was impossible to initially correlate the clinical symptoms with the LORICRIN genotype. Since we have shown that the c.10C > T (p.Gln4Ter) variant was of paternal origin, the father was clinically evaluated. There was no history of skin involvement, but dystrophic nails and massive carries from the age of 20 were reported. Nail involvement in LK is uncommon and also was not described in knockout mice models [1,22,24], although loricrin is expressed in the nail proximal fold [25]. Nevertheless, considering the fact that the father of family 2's history of dystrophic nails was negative, as well as the fact that there were no nail symptoms in the ARCI-affected proband, the nail dystrophy of the father seemed to occur independently. There were also no cases of massive caries among the father's relatives. Interestingly, previous studies have shown that in murine and human aggressive periodontitis, LORICRIN mRNA expression was diminished [26,27]. Therefore, though no skin changes were noted, an open question remains as to whether the presence of the heterozygous PTC variant confers susceptibility to caries.
In conclusion, our results broaden the knowledge about LORICRIN gene variants and their phenotypic significance and give insight into the molecular pathology of loricrin keratoderma lesions.

Skin Biopsy
The transcriptome analysis was performed using RNA isolated from lesional epidermis. The 3 mm skin biopsy from the lesion located on the upper tibia was taken from the LK patient and from the same location of a healthy age-matched male. The biopsies were immediately frozen and kept at −80 • C. The epidermis was mechanically detached from the underlying skin layers in a cryotome prior to RNA isolation.

RNA Sequencing
The samples were mechanically homogenized, and RNA was isolated using an RNeasy Micro Kit (Qiagen, Hilden, Germany). The quality and integrity of total RNA were assessed with an Agilent 2100 Bioanalyzer using an RNA 6000 Pico Kit (Agilent Technologies, Ltd. Santa Clara, CA, USA) In total, polyA enriched RNA libraries were prepared using the QuantSeq 3 mRNA-Seq Library Prep Kit according to the manufacturer's protocol (Lexogen GmbH, Vienna, Austria). Briefly, libraries were prepared from 5 ng of total RNA. The first step in the procedure was a first-strand cDNA synthesis using reverse transcription with oligodT primers. Then, all remaining RNA was removed to what was essential for an efficient second-strand synthesis. The second-strand synthesis was performed to generate double-stranded cDNA (dsDNA). It was initiated by a random primer containing an Illumina-compatible linker sequence. The obtained cDNA was purified using magnetic beads to remove all reaction components. cDNA libraries were amplified by PCR using starters provided by a producer. The library evaluation was completed with an Agilent 2100 Bioanalyzer using the Agilent DNA High Sensitivity chip (Agilent Technologies, Ltd., Santa Clara, CA, USA). The mean library size was 220 bp. Libraries were quantified using a Quantus fluorometer and QuantiFluor double-stranded DNA System (Promega, Madison, WI, USA). Libraries were run in the rapid run flow cell and were single-end sequenced (75 bp) on HiSeq 1500 (Illumina, San Diego, CA, USA).

Statistical Analysis
The quality of sequencing data was firstly checked with the FastQC program [28]. Then, data were mapped to the reference human genome GRCh38 with a star aligner [29]. The calculation of read counts was performed with the HT seq [30]. All genes with very low expression (below 5 counts) across the examined samples were discarded. Due to no replicates for the differential expression analysis, the edgeR method [31], recommended for such an experimental design, was used. We used the value 0.75 as an approximation of the dispersion parameter based on our previous experience with similar data. The gene ontology was performed using the system PipeR package [32]. As important genes, those with an absolute value of log fold change higher than 3 and an abundance of read measured by log counts per million higher than one were chosen. All statistical analyses were carried out using R software v. 4.2.3 [33].  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to ethical reasons.